Magnetization jump in the XXZ chain with next-nearest-neighbor exchange 
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We study the dependence of the magnetization M with magnetic field B at zero temperature 
in the spin-1/2 XXZ chain with nearest-neighbor (NN) J\ and next-NN J2 exchange interac- 
tions, with anisotropies Ai and A2 respectively. The region of parameters for which a jump in 
M(B) exists is studied using numerical diagonalization, and analytical results for two magnons on 
a ferromagnetic background in the thermodynamic limit. We find a line in the parameter space 
(J2/J1, A1/J1, A2/J1) (determined by two simple equations) at which the ground state is highly 
degenerate. M(B) has a jump near this line, and at or near the isotropic case with ferromagnetic Ji 
and antiferromagnetic J2, with | J2/J1I ~ 1/4. These results are relevant for some systems containing 
CuO chains with edge-sharing CuCU units. 
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I. INTRODUCTION 

In recent years quantum spin chains and ladders have 
been a subject of great experimental and theoretical in- 
terest. Quasi one-dimensional magnetic systems have 
been identified and studied experimentally fl|-|4|, and 
new theoretical studies of the spin-1/2 XXZ chain with 
nearest-neighbor (NN) and next-NN exchange coupling 
(which is equivalent to a two-leg zig-zag ladder) were pre- 
sented @-Q. The Hamiltonian is: 



H = ^[Jx(SfSf +l + SfSf +1 + AxSfS^i) 
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One of the properties of this model is that the magnetiza- 
tion as a function of applied magnetic field M(B) displays 
a jump for certain parameters pUO]. This phenomenon 
called mctamagnetic transition was observed for example 
in TmSe |l|, Pe x Mn 1 _ x Ti0 3 @, Tb!_ x S C:c Mn 2 @, 
and the quasi one-dimensional compound Ba3Cu204Cl2 

PI- 

Previous studies of metamagnetism in the model were 
restricted to Ai = A2 = A Ppd], The results show that 
a jump in M(B) is not possible if A > A c . Following 
the methods explained in section III, we have determined 
A c = (-5 + VT7)/4 = -0.22 |9j. In principle, the re- 
quirement of an opposite sign for the z component of the 
exchange than for the other two seems unrealistic. How- 
ever, for the Hamiltonian Eq. (|l|), a negative Ai with 
positive J\ is equivalent to a positive Ai with negative 
J\ , since a rotation of every second spin in 7r around the 
z axis changes the sign of the x and y components of 
J\. Thus, Ai < A c does not necessarily mean a large 
anisotropy of Ai. Instead, a negative A2 < A c is a very 
large anisotropy of J2 and seems difficult to find in real 
systems. 



In this work we extend the previous searches of meta- 
magnetic transitions to arbitrary values of Ai and A2, 
concentrating our study in |Ai| < 1 and < A2 < 1, 
and particularly near the isotropic case of ferromag- 
netic J\ and antiferromagnetic J2. Several systems con- 
taining CuO chains with edge-sharing CUO4 units, like 
La6CasCu2404i, Li2Cu02, and Ca2Y2Cu50io, are ex- 
pected to lie near this limit |2f|[ . In these systems, the 
Cu-O-Cu angle 9 of the CuO chains is near 90° and there- 
fore, the usual antiferromagnetic NN exchange is largely 
frustrated. As a consequence of virtual processes in which 
the Hund rules exchange integral at the O atoms play a 
role, J\ becomes small and ferromagnetic. We estimate 
the critical field B c at which a jump or an abrupt increase 
in M(B) should exist in these compounds. 

The outline of the paper is as follows. In section II we 
explain how we determine the regions of parameters for 
which a metamagnetic transition is expected. Section III 
contains analytical results for the onset of bound states 
of two magnons in a ferromagnetic background. These 
results and numerical ones in chains of 20 sites are used 
in section IV to present phase diagrams in which the 
boundaries of the metamagnetic regions are shown. This 
section includes numerical results for the ground state 
energy per site as a function of magnetic field E(M), 
curves of M(B), and critical field B c as a function of 
a = J2/J1 for A 2 = — Ai = 1. Section V contains a 
more detailed numerical study of this isotropic case for 
a slightly larger than 1/4. Section VI is a summary and 
discussion. 



II. CONDITIONS FOR THE EXISTENCE OF A 
JUMP IN M(B) 

For the sake of clarity, we anticipate some of our nu- 
merical results for the dependence of the ground state en- 
ergy per site E as a function of magnetization M — S z / L, 



where S z = Y^ Sf is the z component of the total spin 
and L is the number of sites. They are shown in Fig. 1. A 
large portion of the curve E(M) has negative curvature. 
These points are actually not accessible thermodynami- 
cally. The dashed straight line is tangent to E(M) at the 
two points Mi = 0.108 and M 2 — 0.5 (the Maxwell con- 
struction). For all M in the interval (Mi, M2), it is ener- 
getically more favorable for the system to phase separate 
into a fraction x = (M — Mi)/(M2 — Mi) with magneti- 
zation M 2 and a fraction 1 — x with magnetization Mi . 
The energy of the mixture is represented by the dashed 
line. If a magnetic field B is applied to the system, the 
equilibrium magnetization is determined minimizing the 
Helmholz transform G = E — MB, what leads to the con- 
dition B = dE/dM . From inspection of Fig. 1, we see 
that M(B) increases with B, until it reaches the critical 
value B c = \E(M 2 ) - E(M 1 )]/(M 2 - Mi) (= 0.00277 Ji 
for the case of Fig. 1). At B = B c , M(B) suddenly 
jumps from Mi to M 2 . If M 2 is smaller than the satura- 
tion magnetization, M increases further for B > B c , but 
we have not found this situation in the present model. 
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A-f do not show a significant size dependence and it is 
accurately determined in chains of 18 sites. Another 
critical value A a (a) was obtained from the condition 
d 2 E/dM 2 \ M=1/2 = 0. For A > A a the curvature 
d 2 E/dM 2 is positive for all M. The curvature at M 2 
was calculated numerically using: 

d 2 E/dM 2 \ M=1/2 = lim L 2 [E{l/2) + E{l/2 - 2/L) 

L— »oc 



-2E(l/2-l/L)}, 



(2) 



in a periodic chain. In contrast to the case of A-^, and 
even taking L = 50, the result has some finite-size effects 
|9| . This is in part due to the fact (explained in the next 
section) that the ground state near M = 1/2 becomes 
incommensurate for sufficiently large a, and these wave 
vectors cannot be represented in small chains if periodic 
boundary conditions are used J12| . From the numerical 
solution of the problem of two spin excitations on the 
ferromagnetic state for L — ► 00, more accurate values of 
A a (a) were obtained recently for a < 1/2 JL1J]. In the 
region of the (a, A) plane where A^(a) < A < A a (a) a 
metamagnetic transition occurs in the model pjl l| . 

Relaxing the condition Ai = A 2 = A, and keeping for 
example A 2 fixed, limiting values of Ai (A{(a, A 2 ) and 
Af(a, A 2 )) can be defined in the same way and we ex- 
pect in principle that magnetization jumps will be found 
if A{(a, A 2 ) < Ai < Af(a, A 2 ). The same is true in- 
terchanging Ai and A 2 . Fortunately, as described in 
the next section, the A" are given by analytical expres- 
sions, which are very simple near a = 1/4. The A{ are 
determined by numerical diagonalization of chains with 
L = 20 and the results are given in section IV. 



III. THE TWO-MAGNON PROBLEM 



FIG. 1. Energy per site as a function of total z component 
of spin per site for a chain of 20 sites (solid circles). The 
full line is a polynomial fit. Dashed line and diamonds corre- 
spond to the Maxwell construction. Parameters are J\ = 1, 
a = J2/J1 = 0.25, Ai = -0.9 and A 2 = 0.6. 

From the above description, one can see that for 
a metamagnetic transition to occur at very low tem- 
peratures, E(M) should satisfy two conditions: 1) 
d 2 E/dM 2 < in a finite interval of values of M. 2) 
E(M 2 ) > E(Mi), where Mi < M 2 are determined by 
the Maxwell construction. From the general behavior of 
E(M) for the case Ai = A 2 = A , Gerhardt et al. §] have 
found that when metamagnetism exists, M 2 = 1/2 and 
the condition 2 ceases to be satisfied when Mi = 0. More 
precisely, from their finite-size results for E{M, a, A), 
with a = J 2 /Ji, they obtained a critical value of A 
(A* (a)) from the equation E(0,a,Af) = E(l/2,a,A f ). 
For A < A-f the system is a fully spin polarized fer- 
romagnet even at B = 0. Fortunately, the results for 



In this section, we derive analytical expressions for the 
upper boundaries of the expected metamagnetic region 
(Ai(a,A 2 ) or A 2 (a, Ai)), looking for bound states of 
two spin flips on a ferromagnetic background in an in- 
finite chain. This is a two-body problem which due to 
translational invariance can be mapped into a single- 
body one. Similar problems were solved for example to 
determine exactly metal-insulator boundaries in general- 
ized Hubbard models |2M. For Ai = A 2 = A, this prob- 
lem has been solved numerically for L ~ 150 by Cabra 
et al. S and for L — > 00 and a < 1/2 by Hirata [O. 
Simple analytical expressions are given in Ref. |1S|| . 

We use a Jordan- Wigner transformation S^ = 

c]exp(inY, l<:j ni), Sj~ = (#+)*, Sj = n,- - 1, with 



n j = Ca cj , to express the spin operators in terms of spin- 
less fermions. Calling a = J2/J1, setting J\ — 1 as the 
unit of energy, and substracting the constant LE{\/2) 
(with £(1/2) = (Ai + aA 2 )/4), the Hamiltonian Eq. @ 
takes the form: 



!/ t 
i)ni + -{cl +1 d + < 



H = 2J[(-Ai - aA 2 )ni + ^(c\ +1 a + ac J i+2 a + H.c. 

i 

+A 1 n i n i+1 + aA 2 n i n i+2 -a(4 +2 n i+ ic i + R.c.)]. (3) lill( , a ,. hl ( , )S , y an(1 llsiu; , [22 



from q = tt for a < 1/4 to g = ± arccos(— l/(4a)) for 
a > 1/4. The integrals can be solved decomposing the 
integrand into a sum of expressions with denominators 



After Fourier transform Cj = \/l/L ^2 e lq ^c q , the model 
can be written as: 

H = '^^e q c) q c q + — ^ ^2 [(Ai + 2a cos X) sing sing' 



-aA 2 sm2qsm2q'y & _ ,c'K +q ,CK +q CK^ q , 



K q,q'>0 

■,nJ c t 

-q' 



(4) 



with 



- Ai — aA 2 + cos q + a cos 2q. 



The two-magnon eigenstates of total momentum K can 
be written in the form \ip(K)) = Y^, q> o A c *L_ g c4 + j0). 

Replacing this into the Schrodinger equation H\ip) — 
\\ip), one obtains: 

-4[(Ai + 2acosK)S 1 (K) sing + aA 2 S 2 (K) sin2g] 



£k. 



^+9 



where 



1 



S n (K) = j^A q sinng. 



(5) 



(6) 



Replacing Eq. (||) into Eq. (JsJ) leads to two homoge- 
neous equations for the two unknown sums Si and 5*2. 
For fixed values of the parameters and K, the condition 
of vanishing secular determinant determines the allowed 
eigenvalues A. Since we are looking for the condition on 
the parameters for the onset of a bound state (when the 
right hand side of Eq. (En becomes negative), we set A 
slightly below two times the minimum one-magnon en- 
ergy: A = 2 min(eg) — 77, where 77 is a positive infinitesimal 
energy. With this value of A, the condition of vanishing 
secular determinant in the thermodynamic limit L — > oo 
takes the form: 

4(Ai + 2acos#)[Io + 16aA 2 (I I 2 - if)} 

+ 16aA 2 I 2 + 1 = 0, (7) 

where the integrals I n (a, K) (n = 0, 1, 2) are: 



In = 7T- 



dq sin q cos™ q 



27r 7o 4a cos K cos 2 q + 2 cos -y cos q + c(a, K) + r\ ' 

(8) 

with r\ — > and 

c(a,K) = 2(l-a-acosK), a < 1/4 

c(a, K) = 2a H 2a cos K, a> 1/4. (9) 

4a 

The change of the expression for c(a,-ftT) at a = 1/4 is 
due to the change in the wave vector which minimizes e q 



1 

2^ 



dq sin < 



b cos g 




(10) 



For each value oi a, K and Ai, Eq. (Q) is a linear equa- 
tion in A2. The searched upper boundary A§(a, Ai) is 
determined choosing the value of K which leads to the 
highest root of Eq. ([!]). For A 2 < A 2 (a, Ai), the curva- 
ture Eq. (0) is negative. The same is true interchanging 
Ai and A 2 . If Ai = A 2 = A is taken @|I|], the highest 
root of the quadratic Eq. (pi) determines A a (a), since 
there is at least one bound state for A < A a (a). 

For a < 1/4, the total wave vector K which first leads 
to a bound state is K = 0. In the sector of two particles, 
this is also the wave vector of the ground state of the 
non-interacting part of the Hamiltonian in the fermionic 
representation (Eq. (||) or (|J)). Using Eqs. (||), @ and 
(p0|), we obtain after some algebra: 



la 



r-1 l-(l-2a)r 

h — 



8a 



16a 2 



a — 1 



1 



-1 /( 



16a 2 ' \2a 
with r = (l-4a)~ 1/2 , (a < 1/4). (11) 

Eqs. (0) and @ define Af(a,Ai) and A?(a,A 2 ) for 
a < 1/4. 

For a > 1/4, the ground state of the non-interacting 
fermionic Hamiltonian for two particles is degenerate 
with total wave vector K = or K = ±Ki with Ki — 
2arccos(— l/(4a)). For Ai = A 2 , we have found that 
the wave vector of the ground state of the two-magnon 
problem for parameters near the onset of a bound state is 
K = for a < 1/2. At a = 1/2 it jumps to ±K t . From 
numerical investigations in finite systems (with L ~ 20), 
using twisted boundary conditions to allow all possible 
wave vectors Jl2| , we find that for values of a > a w ~ 1 , 
the ground state wave vector deviates continuously from 
±i-Q. Comparison of the numerical results for A a (a) 
(with L ~ 50) of Ref. J9| with our analytical ones assum- 
ing K = ±Ki, gives a w = 0.77. In the region a > a w , 
it seems not possible to find analytical results for A a (a). 
However, taking K = Ki leads to a reasonable lower 
bound, since K is not too different from Ki (in any case, 
as we shall see in section V, Eq. (g) ceases to be valid for 
ot, > a w ). In the general case Ai ^ A 2 , the maximum be- 
tween the results for A„ assuming K — or K — Ki, also 
gives a lower bound for A 2 (a, Ai) and Af(a, A 2 ). We 
restrict the calculation to these values of K. For K = 
and a > 1/4, the integrals Eq. (||) diverge for 77 — > 0, 
and Eq. ( |l0| ) has complex coefficients for finite 77. The 
equation for the A„ is obtained after a careful limiting 
procedure of Eq. (0) with adequate choice of the branch 



of the root in Eq. (|l0[) . Physically this means simply 
to consider values of A n slightly smaller than the critical 
ones, in such a way that they lead to a finite binding en- 
ergy r], and then take the limit 77 — > 0. The final results 
turns out to be very simple: 



AiA 2 + Ai 
2a 



A 2(1 + ^T 



+ 1 = 0, 
(a > 1/4, K 



0). (12) 



As an example, this equation is satisfied for a = 1/4, 
Ai = -0.9, A 2 = 2/3. Lowering A 2 a little bit, E(M) 
displays the behavior required for a jump in M(B) (see 
Fig. 1). 

Finally, for a > 1/4 and K — Ki, the integrals take 
the form: 



h = 



1 



1- 



1 



8aB V Aay/A 
l/l 1 2a 



,h = 



1 



1 



1 



32(aB) 2 \2a y/J 



A 



16a \ A B 2aAB 3 



, with 



A= 1--L 5 = 1- -L ( a > i/4,K = Ki). (13) 
16a z 8a z 

It remains to establish when A" are determined by Eq. 
(H), or Eq. (|l|) and Eq. (@) with cosK = l/(8a 2 ) - 1. 
Taking for instance Ai fixed, equating the values of A 2 
obtained from both expressions, we obtain: 

= -1 + 2A?, (-1 < Ai <0). (14) 



4Ai 



For smaller values of a, AS^a, Ai) is determined by Eq. 
(O). Similarly, for a given value of A 2 , Eq. (U2) gives 
A?(a,A 2 ) if a < J(l + A 2 }/2]- 1 / 2 /4, while for larger 
values of a, Eqs. (|7) and ( |l3| ) should be used. 

Eq. ( \\A\) defines a particular line in the parameter 
space (a, Ai, A 2 ), which as will be seen in the next sec- 
tion, plays an important role in our study. From the ana- 
lytical results of this section, we see that at this line, not 
only there is a degeneracy in the lowest lying eigenstates 
in the two-magnon sector (for M = 1/2 — 2/L with L — > 
00) between the wave vectors K — and K — ±Ki, but 
also lim^oo L(E(l/2 - l/L) - £(1/2)) = e arccos(Al ) = 
and limL^oo 1(^(1/2 - 2/L) - £(1/2)) = 0. In other 
words E(M) is independent of M for M ~ 1/2 in the 
thermodynamic limit. 



IV. PHASE DIAGRAM FOR METAMAGNETIC 
TRANSITIONS 

In this section, we delimit the region of parameters 
inside which a jump in the function M(B) is expected. 
Specifically for each a and Ai, we determine by numeri- 
cal diagonalization of chains of 20 sites the lower bound- 
ary A2 of the metamagnetic region (from _E(l/2) = 
E(0)), while the upper boundary A2 (or a lower bound of 



it for large a) is taken from the expressions of the previ- 
ous section. The same is done interchanging Ai and A 2 . 
This study is complemented by numerical calculations of 
the ground-state energy as a function of magnetization 
E(M) (from which M(B) can be derived), and the crit- 
ical field B c for several parameters inside the region of 
interest. 




FIG. 2. Boundaries of the region of expected metamag- 
netic transitions in the (a,A2) plane, for different values of 
Ai. This region is shown dashed for Ai = —0.7. The symbols 
(joined by thin lines) indicate the lower boundary Aj calcu- 
lated numerically in chains of 20 sites. The upper boundary 
A§ (thick lines) is taken from the analytical results of section 
III. 

In Fig. 2 we show A2 (discrete points) and A 2 (thick 
curve above it) as a function of a for several values of 
Ai. The region of interest A 2 < A 2 < Aj is shown 
by vertical dashed lines for Ai — —0.7. Clearly the up- 
per boundary A2 displays two kinks as a result of the 
change in the ground state wave vector for one magnon 
at a = 1/4, or two magnons at a = — l/(4Ai), as de- 
scribed in the previous section. A 2 is given by Eqs. ( pd| ) 
and (0) with K = for a < 1/4, by Eq. (|lg) between 
both kinks, and (as a lower bound) by Eqs. (O) and 
(0) with cos if = l/(8a 2 ) - 1 for a > -l/(4Ai) respec- 
tively. It is remarkable that at the second kink (deter- 
mined by Eqs. (H)), where the wave vector of the two- 
magnon ground state jumps from K = to K — ±Ki 
and dE/dM\ M=1/2 = d 2 E/dM 2 \ M=1/2 = for L -> 00 
according to the results of the previous section, also 
E(l/2) — E(0) (at least within the accuracy of the Lanc- 
zos diagonalization an independently of system size) . Ac- 
tually E(M) is quite flat at this point. Specifically, wc 
find by numerical diagonalization that if and only if the 
energy for each value of S z is minimized with respect 
to the optimum twisted boundary conditions ]12), E(M) 
becomes independent of M at this point, within our pre- 
cision of 10 -9 . All wave vectors at the minimum become 
incommensurate except for S z — and S z = L/2. For 
S z = L/2 — 2, we find that the minimum is at wave vec- 



tor K = Ki, while the energy at K — is of the order 
of 1CP 5 above the minimum energy for L = 20, due to 
finite size effects. 

Particular cases of degeneracy at the line determined 
by Eqs. (Oh are already known. The degeneracy at the 
point a = 1/4, Ai = — 1, A 2 — 1 was studied by Hamada 
et ai, who also have shown that the ground state in 
the sector of total spin S = is a resonance- valence- 
bond state involving singlet pairs at all distances p3[ . 
The degeneracy at the Majumdar- Ghosh point a = 1/2, 
A x = A 2 = —1/2, was discussed by Gerhardt et al. 0. 
For Ai — > 0, this point moves towards two decoupled fer- 
romagnetic Heisenberg chains (see Eqs. (14)), for which 
E(M) is constant. 

This special point where A^ and Aj coincide separates 
the region of expected metamagnetism in two zones. The 
one at the left of the point shrinks and moves towards 
A 2 = 1 as Ai decreases approaching the isotropic limit 
Ai = — 1, where it disappears. The right zone also dis- 
places towards A2 = 1, but increases as Ai moves to 
— 1, suggesting that a jump in M(B) is also possible in 
the isotropic limit A 2 = — Ai = 1, if a > 1/4. Fig. 3 
shows the same two metamagnetic zones in the (a,Ai) 
plane for several values of A 2 . The same tendencies as 
before are observed as the isotropic limit A 2 = — Ai = 1 
is approached. 



-0.80 





FIG. 3. Boundaries of the region of expected metamag- 
netic transitions in the (a, Ai) plane, for different values of 
A2. This region is shown dashed for A2 = 0.5. The symbols 
(joined by thin lines) indicate the lower boundary A{ calcu- 
lated numerically in chains of 20 sites. The upper boundary 
Af (thick lines) is taken from the analytical results of section 
III. 

In the zone of lower a, the existence of a jump in M{B) 
is confirmed by numerical calculation of E(M). An ex- 
ample was shown in Fig. 1. Instead, inside the right zone 
of expected metamagnetism (for a larger than the point 
of degeneracy at which A2 and AjJ meet), the situation 
is not so clear. As it is clear in Fig. 4 for a — 1/2, the 



FIG. 4. Energy per site as a function of total spin per site 
for a chain of 20 sites with J\ — 1, Ai = — 1, A2 = 1, and 
several values of a — J2/J1. Full lines are polynomial fits (see 
text). Dashed line and diamonds correspond to the Maxwell 
construction. 

energy per site as a function of the z component of the 
total spin S z = Yli S? — ML, shows a significant even- 
odd effect for small chains. This effect has been reduced 
minimizing E{M) with respect to the optimum twisted 
boundary conditions to allow for incommensurate wave 



vectors. |1J] Most of the resulting wave vectors are in- 
commensurate, particularly for odd S z and low a > 1/4. 
In spite of this procedure, for a > 0.4, all E(M) for 
odd S z seem to be shifted to higher energies in compar- 
ison with the corresponding values for even S z . If this 
tendency persists in the thermodynamic limit (keeping L 
even) states with odd S z would not be accessible thermo- 
dynamically, and the calculation of the previous section 
(based on Eq. (||)) would be irrelevant. In any case E(M) 
is very flat for 0.3 < M < 0.5 (see Fig. 4) and if the cur- 
vature at M = 1/2 is positive, there would be a steep 
increase in M{B) for M > 0.3, which is probably hard to 
distinguish experimentally from a jump. The determina- 
tion of the critical value of a for which d 2 E/dM 2 \ M =i/ 2 
changes sign, is a difficult task which is postponed to the 
next section. 



Li 2 Cu0 2 , and Ca 2 Y 2 Cu50io, 
and 65 Tesla respectively. 



20] we obtain B c = 14, 40 
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FIG. 5. Total spin per site as a function of magnetic field 
for J\ — 1, Ai = — 1, A2 = 1, and two values of a. 

The continuous curves in Fig. 4 correspond to fits in 
the numerical data using a polynomial of even powers 
of M with about half as many parameters as points to 
be fitted (to average the even-odd effect). These con- 
tinuous curves allowed us to perform the Maxwell con- 
struction analytically (dashed lines and diamonds in Fig. 
4), and to calculate B(M) from B = dE/dM (see sec- 
tion II). The magnetization curve B(M) is shown in Fig. 
5 for two values of a, which are near those estimated 
for La6CasCu 2 404i, and Ca 2 Y 2 CusOio respectively. Q 
While the details of the curve and the value of the critical 
field at the jump B c (if it exists) depend on the fitting 
procedure, the fact that there is an abrupt increase from 
~ 60% to 100% of maximum magnetization with a very 
small variation of magnetic field is a genuine feature of 
the system. To estimate the variation of B c with a for 
the isotropic model A 2 = — Ai = 1, we have calculated 
the average slope oiE(M) (minimized with respect to the 
optimum twisted boundary conditions) between M = 0.3 
and M = 0.5. The result is shown in Fig. 6. Assuming 
.9 = 2 for the gyromagnetic factor of the Cu ions, and tak- 
ing the values of J\ and a estimated for LagCasCu^CUi, 



FIG. 6. Critical magnetic field as a function of a for J\ = 1, 
Ai=-1, A 2 = l. 



V. MAGNETIZATION JUMP IN THE 
ISOTROPIC CASE 

The results of the previous section for A 2 = — Ai = 1, 
show that there is an abrupt increase in M(B) for B ~ 
B c (a) near M = 1/2, particularly for 1/4 < a < 1/2. 
However, they are not enough to establish the existence 
of a true jump. The calculation of d 2 E/dM 2 \ M =i/ 2 be- 
comes complicated by the fact that, independently of sys- 
tem size (at least for L < 40 and near M = 1/2), only 
states with total spin S = L/2 — In min with /, n m ; n inte- 
gers, but n m i n > 1, can be of thermodynamic relevance. 
In other words, the curve E(M) near M = 1/2 looks like 
a straight line of slope B c plus a periodic function with 
period n m i„/ L (see for example Fig. 4 for a = 0.5, where 
n min = 2; and Fig. 7 for a = 0.3, where n min = 3). If 
d 2 E/dM 2 \ M= i/ 2 > 0, this means (as stated in Ref. 8) 
that reducing the magnetic field from values high enough 
to ensure saturation of the magnetization, the spins flip 
in groups of n m i n - The physical reason of this behavior 
is not completely clear. It seems that spin flips tend to 
bind in groups of n nlm . 

As a consequence, Eq. (||) which could be calculated 
analytically becomes invalid and should be replaced by:: 



d 2 E/dM 2 \ M=1/2 = lim -J— [E{l/2)+E(l/2-2n min /L) 

L ^°° n min 

-2E(l/2-n min /L)}. (15) 



While n n 



2 for a > 0.4, 



increases with de- 



creasing a, making larger system sizes necessary for an 
accurate evaluation of Eq. (fla). In addition, this equa- 
tion cannot be used for values of a for which n mul is 
not well defined. For example a ~ 0.35 (see Fig. 4) 
seems to correspond to a transition from n m i n = 2 to 



Hmin = 3 with lowering a. Taking into account these dif- 
ficulties, we have chosen three values of a in the range 
1/4 < a < 0.3, with apparently well defined values of 
«min (3, 4 and 5 for a — 0.3, 0.28 and 0.26 respectively, 
see Fig. 7), and have calculated Eq. (H3) as a function 
of system size. For a = 0.3, for which the energy with 
six (2n m in) spin flips should be calculated, the largest 
system size considered was L = 40. This was reduced to 
L = 28 and L = 24 for a = 0.28 and 0.26 respectively, 
due to the increase in 2n m i n - We have taken L multiple 
of four to avoid frustration of the next-nearest-neighbor 
antifcrromagnetic interaction J 2 . 
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the three values of a are shown in Fig. 7. The oscil- 
lations with period n m i n /L are clearly seen. We also 
show a comparison with the result of periodic bound- 
ary conditions. As before, the minimization with respect 
to twisted boundary conditions reduces the magnitude 
of these oscillations and the finite size effects, partic- 
ularly for a near 1/4, small M and smaller L. With 
increasing L, the points tend to lie nearer to the dot- 
dashed line. As mentioned before 12], for any system 
size, minimization with respect to twisted boundary con- 
ditions reproduces the exact energy and wave vector in 
the one magnon sector: for a > 1/4 and M = 1/2 — 1/L, 
q = ±arccos(— l/(4a)) and 

E(l/2-l/L) = (Ai+oA 2 )Ji/4-[Ai+a(l+A a )+l/(8a)]/£. 

However, for a — 0.3, evaluation of Eq. ( |l5| ) is not af- 
fected by the use of periodic boundary conditions, since 
they minimize the energy for three and six spin flips. 
Then, our results for a = 0.3 are consistent with those of 
Cabra et al., who using periodic boundary conditions ob- 
tain that M(B) is very steep near M = 1/2, but without 
a jump for a = 1/3 (Fig. 4 of Ref. g). 
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FIG. 8. Relative 

curvature C = d 2 E/dM 2 \ M=1/2 /(E(l/2) - £7(0)) calculated 
using Eq. (h_5[) as a function of the inverse of the system size 
for Ji — A2 = — Ai = l, and several values of a. 



FIG. 7. Energy per site as a function of magnetization for 
Ji = A2 = — Ai = 1, and several values of a = J2/J1 
and system sizes: solid circles L = 20, empty diamonds 
L = 24, and empty circles L — 28. Crosses correspond to 
L — 20 using periodic (instead of twisted) boundary condi- 
tions. The dot-dashed line joins the points with M — 1/2 and 
M = 1/2 - 2n min /L for L = 20. 

The results of E(M) for two different system sizes and 



In Fig. 8, we represent the dimcnsionlcss relative cur- 
vature: 



C 



d 2 E/dM\ 1=1/2 
£(1/2) - £7(0) ' 



where the numerator is evaluated using Eq. (|15|) with 
different values of L, and in the denominator the result 
for L = 20 is taken. Except for the case of a = 0.28, in 
which there seems to be an oscillating behavior of C with 
L, the variation of C with L is rather smooth. For L = 24 
(1/L ~ 0.042) the values of C are shifted upwards with 



respect to a smooth curve which fits the rest of the points. 
For a = 0.26, different extrapolations of the data to the 
thermodynamic limit give \C\ < 0.1. Since this value 
is of the order of the uncertainty in the extrapolations, 
a definite conclusion regarding the sign of C cannot be 
drawn. Similary, for a = 0.28, due to the oscillating 
behavior of the data, no reliable extrapolation can be 
made. Instead, for a = 0.26, although we have only two 
points available, it is remakkable that the curvature is 
practically constant (C = —1.617) suggesting a negative 
value in the thermodynamic limit. Taking into account 
that for a — 0.28 and 0.26, the curvature for L = 24 
seems shifted to higher values in comparison with the 
rest of the curve, the sign of C for all finite L, and the 
behavior of E(M) displayed in Fig. 7 for each value of a, 
we believe that C changes sign for a c slightly below 0.28. 
For 0.25 < a < a Cl we expect a true jump in M(B). 



VI. SUMMARY AND DISCUSSION 

We have investigated by analytical and numerical 
methods, the regions of parameters Ai, A2 and a — 
J2/J1 for which a jump in the magnetization as a func- 
tion of magnetic field M(B), of the spin-1/2 XX Z chain 
with next-nearest-neighbor exchange (Eq. (Qj)) is ex- 
pected. The numerical results are restricted to the region 
I Ai| < 1 and < A 2 < l.We were particularly interested 
in parameters near the isotropic limit of ferromagnetic J\ 
and antiferromagnetic J 2 (ot > 0, A 2 = — Aj = 1), which 
are relevant to some systems containing CuO chains with 
edge-sharing CuC>4 units. |2C| | 

One of the necessary conditions for a jump in M(B) 
is that the ground-state energy per site as a function 
of magnetization E(M) should have zero or negative 
curvature in a finite interval. In absence of even-odd 
effects or spin flips in groups as discussed in sections 
IV and V respectively, the curvature at maximum M 
(d 2 E / dM 2 \ M=l / 2 ) can be calculated analytically in the 
thermodynamic limit from the energy of the states with 
one and two magnons. The resulting points of zero curva- 
ture define a surface S\ in the parameter space a, Ai, A 2 , 
which turns out to be a frontier for the occurrence of 
metamagnetic transitions. The rest of the boundary of 
the region of expected metamagnetism was constructed 
from the condition £'(0) = £(1/2), where E(M) was cal- 
culated numerically in chains of 20 sites. We call this 
surface S 2 . 

For — 1 > Ai > 0, the line of Si given by: 

4qAi + 1 = 0, 1 + A 2 - 2A? = 

on which the wave vector of the two-magnon ground state 
in the thermodynamic limit changes from K = to K — 
±2 arccos(— l/(4a)) coincides with a line of S 2 in which 
E(AI) is independent of AI. This noticeable property 
allows to split the regions of expected metamagnetism in 
two, depending if a is moved to lower or higher values 



from this line (see Figs. 2, 3). Inside the first region, 
the existence of a jump in M(B) is confirmed calculating 
E(M) numerically, for all possible M in a finite chain. 
At least for |Aj| < 1, this region disappears if Ai = —1 
or A 2 = 1. However, particularly for a = 1/4, there 
are values of Ai and A 2 near the isotropic case A 2 = 
— Ai = 1 for which a metamagnetic transition exists. 
The intersection of Si with the plane a = 1/4 gives the 
simple equation: 

2AiA 2 + 2Ai + 3A 2 + 1 = 0. 

Displacing one of the Aj slightly to lower values from this 
line (without reaching S 2 ) is enough to have a jump in 
M(B) (see Fig. 1). 

In the second region of expected metamagnetism 
(larger values of a), the analytical calculation of Si be- 
comes invalid due to the tendency of the system to de- 
crease the magnetization from saturation in more than 
one spin flip. Numerical calculations of the curvature 
d 2 E/dM 2 in the isotropic case A 2 = — Ai = 1 near 
AI = 1/2 using states with adequately chosen total spin 
and twisted boundary conditions (which allow for incom- 
mensurate wave vectors |l2|| and are crucial here), suggest 
that it remains negative for a < a c with a c ~ 0.28 or 
slightly less. A jump in M(B) exists for 1/4 < a < ct c . 
In any case, even for a > a c , our results seem enough to 
show that AI(B) should have an abrupt increase (if not a 
true jump) from nearly 60% to 100% of the magnetization 
of saturation at a critical field B c . As an example, taking 
the parameters of Ref. pp] , we estimate B c =14 and 40 
Tesla for La6Ca 8 Cu 2 404i and Li 2 Cu0 2 respectively. 
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